arXiv: 1505.07516vl [astro-ph.GA] 28 May 2015 


Draft version May 29, 2015 

Preprint typeset using style emulateapj v. 5/2/11 


GRAVITATIONAL ENCOUNTERS AND THE EVOLUTION OF GALACTIC NUCLEI. 1. METHOD 

David Merritt 

School of Physics and Astronomy and Center for Computational Relativity and Gravitation, Rochester Institute of Technology, 

Rochester, NY 14623 
Draft version May 29, 2015 

ABSTRACT 

An algorithm is described for evolving the phase-space density of stars or compact objects around 
a massive black hole at the center of a galaxy. The technique is based on numerical integration of the 
Fokker-Planck equation in energy-angular momentum space, f{E, L, t), and includes, for the first time, 
diffusion coefficients that describe the effects of both random and correlated encounters (“resonant 
relaxation”), as well as energy loss due to emission of gravitational waves. Destruction or loss of stars 
into the black hole are treated by means of a detailed boundary-layer analysis. Performance of the 
algorithm is illustrated by calculating two-dimensional, time-dependent and steady-state distribution 
functions and their corresponding loss rates. 

Subject headings: galaxies: evolution — galaxies: kinematics and dynamics — galaxies: nuclei 


1. INTRODUCTION 

The distribution of stars around a massive black hole is a well-studied problem. Many solutions are consistent with 
Jeans’s theorem, but after a sufficiently long time, gravitational encounters between stars will drive their distribution 
toward a predictable form. In the case of stars of a single mass, random encounters result in a steady-state distribution 
that was hrst correctly described by Bahcall and Wolf (1976): the density of stars falls off as n(r) oc and the 

phase-space density is f(E) oc \E\^^^] here r is distance from the black hole and E is the orbital energy per unit mass 
of a star. Bahcall and Wolf’s derivation assumed that the gravitational potential was dominated by the black hole, 
and that the distribution of stars in velocity space was isotropic. Capture of stars by the hole was allowed to occur 
only via diffusion in energy; the steady-state feeding rate was found to be small, in the sense that a small fraction of 
stars, within any radius, are consumed by the hole in one, two-body relaxation time at that radius. For this reason, 
/ ^ is often described as a “zero flux” solution and in fact it can most easily be derived by setting to zero the 

e nergy-space flu x in th e evolution equation for the orbital distribution. 

iFrank fc ReesI (II976I1 argued that feeding of the black hole would be dominated by diffusion in angular momentum, 
no t energy, and that a f ractio n of order unity of sta r s with in radius r would be captured in one relaxation time at 
r. iLightman fc Shapirol (jl977ll and iCohn fc Kiilsrudl (jl978l l found steady-state solutions of the equations describing 
diffusion in both E and L; capture or destruction of stars was assumed to occur when the orbital eccentricity was 
large enough, at a given energy, that the orbit intersected the loss sphere, the region around the hole in which stars 
are consumed or tidally disrupted. The steady-state solutions found by these authors were anisotropic with respect to 
velocity, since the phase space density must be zero near the loss sphere. However the steady-state energy distribution 
was found to be quite similar to that of the scale-free Bahcall-Wolf solution, at least for energies much greater than 
that of a circular orbit at the edge of the capture sphere, and the corresponding configuration-space density was close 
to n ^ 7.-774 outside of this sphere. 

These early treatments were made possible by the availability of analytic expressions for the diffusion coefficients that 
describe encounter-driven changes in E and L. Such expressions can be straightforwardlv d erived in the case of random 
encounters in a homogeneous medium fe.g. iRosenbluth et al.llI957l : ICohn &: Kulsrudlll978ll . Sufficiently near to a mas¬ 
sive black hole - very roughly, at distances less than ^ 10 “^ times the black hole’s gravitational influence radius - the 
assumption of random encounters breaks down, since orbits resemble closed Keplerian ellipses for many periods. In this 
regime, diffusion in angular momentum is dominated by “resonant relaxation” (RRl (|Rauch fc Tremain^ll996fl while 
changes in energy are still controlled by random (“non-resonant,” NR) encounters. Approximate timescales for changes 
in L in the RR regime have long been available, but until recently, usefully accur ate expressions for the diffusion coef¬ 
ficien ts in L were not available. Happily, that situation has now begun to change dHamers. Portegies Zwart fc Merritt! 
I2014fl , and so it is feasible to repeat calculations like those of iCohn fc Kiilsriidl (Il978ll using expressions for the diffusion 
coefficients that are valid much closer to the black hole. 

This paper, which is the first in a series, presents a numerical algorithm for solving the Fokker-Planck equation 
describing the evolution due to gravitational encounters of f{E,L), the two-dimensional ph ase-space density of star s 
orbiting around a massive black hole. The basic numerical approach is the same as that of ICohn fc KulsrudI (I1978I1 . 
although their algorithm has been generalized to include diffusion coefficients that describe both random and correlated 
encounters. Cohn & Kulsrud’s treatment is improved upon in other ways as well; notably by the use of a logarithmic 
grid in angular momentum, which is necessary to accurately treat the behavior of highly eccentric orbits. 
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As far as we are aware, the time-dependent, f{E,L,t) solutions presented here are the first ever published. While 
Cohn & Kulsrud’s algorithm was capable of calculating time-dependent solutions, those authors (perhaps because of 
computer limitations) chose to present only steady-state solutions in their 1978 paper. Time dependence is particularly 
relevant to galactic nuclei since (energy) relax ation times ar e believed to be comparable with galaxy lifetimes, even 
in nuclei as dense as that of the Milky Way ([Merrittl 1201011. This time dependen ce has routinely been ignored in 
calculations of the rate of stellar tidal disruptions ie.g. iMasrorrian fc Tremain^ 1 19981 1. 

In spite of the improvements, the algorithm presented here still contains some basic limitations. The mass of the 
black hole is fixed and its spin is ignored. Stars are assumed to have a single mass, and the contribution of the stars to 
the gravitational potential is ignored. As a consequence, results are limited in their applicability to a region inside the 
black hole’s sphere of gravitational influence, and the possible influence of spatial asymmetries in the stellar potential 
on the behavior of orbits is ignored. These is no “source term” corresponding to star formation and binary stars are 
not allowed. Rotation of the stellar cluster is ignored. Some of these restrictions will be lifted in later papers from 
this series. 


2. ASSUMPTIONS AND BASIC RELATIONS 

Stars - a term used here to refer both to normal stars and to compact objects, e.g. stellar-mass black holes - are 
assumed to have a single mass, m*. The stars are assumed to be close enough to the black hole (SBH) that the 
gravitational potential defining their unperturbed orbits is constant in time and given by 


<i)(r) = —ipir) = — 


GM, 


( 1 ) 


+ = L = \rxv\. 

Following iCohn fc KulsrudI (1197811 we define new variables {£,TZ) as 

y2 t2 

£ = -E = -—+tP{r), T^=j 2 
with Lc{£) the angular momentum of a circular orbit of energy £\ 


with M, the SBH mass, also assumed constant. The contribution to the potential from the stars themselves is ignored, 
hence results are only expected to be accurate inside the SBH gravitational influence sphere. 

Unperturbed orbits respect the two isolating integrals E (energy) and L (angular momentum), both defined per unit 
mass: 

„,2 „,2 r<i\/r 

( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 

( 6 ) 

( 7 ) 


Ll{£) = 


(GM,) 

2 £ 


Hence 0 < 7^ < 1. £ and TZ are related to the semimajor axis a and eccentricity e of the Keplerian orbit via 

GM, 2 


The orbital (Kepler) period is 


and the gravitational radius of the SBH is 


P{a) = 


^ = 1 - 7^. 

2 £ ’ 

27ra^/^ TT GM, 


y/GW, y/2 f3/2 


r„ = « 4.80 X 10 


-5 ( 

\10^ Mq 


mpc. 


Spin of the SBH is ignored. 

Stars are assumed to be lost - either captured, or destroyed - if their (Newtonian) periapsis distance falls below ric, 
i.e. if 

a(l - e) < He. (8) 

For stars not subject to tidal disruption, e.g. compact objects, rje ~ 8rg (|Merritfll2013L §4.6); otherwise Hc > 8rg is 


the tidal disruption radius. Defining 


£ic = 


GM, 

2nc 


the energy of a circular orbit at r = ric, allows the loss condition to be written as 

7^<7^lc(f) = 2^(l-i^), £<£ic. 


£ic 


2£ic 


( 9 ) 


( 10 ) 
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For E fic, equations m and (US imply 

£ £ 

7^lc 2- « 32^; 11 

£ ic c‘^ 

the final expression assumes £\c = c^/16, appropriate for ric = 8rg. 

The number density of stars in phase space, /, is assumed to satisfy Jeans’s theorem at any given time, but the 
dependence of / on £1 and TZ is allowed to vary with time as gravitational encounters cause stars to change their orbits: 

f = f{£,n,t). ( 12 ) 


The configuration-space density is given in terms of / as 


n{r) = ^ 


ftp{r) fU, 

Ll{£) d£ 


f{£,TZ)dn 


/o 


= V2'; 


GM, 


0 \/2 [V^(r) -£]- nLl{£, t)/r'^ 

d£ /(£:,7^,^)d7^ 

7^ Jo 




where 


TZ„ 


^ 2r^ - £] ^ / _ 

^ m£) [gmJ 1 , 


r£ 

gmI 


(13a) 

(13b) 


(14) 


with similar expressions for the velocity dispersions Or, ot in the radial and transverse directions. 

The time dependence of / is assumed to be described by the orbit-averaged Fokker-Planck equation: 

BN r) 1 B 1 B"^ 

^ + 2 ^ (^((^^) )t) - ^ iN{ATZ)t) + (NiiAn) )t) + {N{A£An)t) 

(15) 


where 


N{£,n,t)=ATr‘^Ll{£)P{£,n)f{£,n,t) = V2 tt^ iGM,f f{£,TZ,t) = J{£,TZ)f{£,'R.,t) (16) 

is the distribution of orbital integrals. The quantities in () in equation (fT^ are diffusion coefficients; the subscript “t” 
indicates a time average over the unperturbed orbit. The diffusion coefficients are functions of £, TZ and of the stellar 
distribution itself, i.e. of /; their functional forms are discuss ed in more detail below . 

Equation (1X5]) can be written in flux-conservation form as (jCohn fc Kulsrudlll978ll 


d f d d 

-4>s = Dee-J^ + Deh-^ + Dsf, 
BE BTZ 

—4>n = Dns-J^ + Dnn-^ + D-jif 
BE BTZ 


with “flux coefficients” 


D, = -(AO, - + ~(ACAK). • 

Dn = -{ATZ)t - ^{A£ATZ)t + \^{A£ATZ)t + , 

De£ = ^{{A£f)t ,Den = D^s = = ^((ATe)')* . 

The rate of loss of stars past the loss-cone boundary: 

TZ = TZici£), £^min A £ A Emax 

is given in terms of the fluxes as 


N = 


d£ 


^min 




d£ 


J-R-ioiS) 

J{£) {£)]d£- 


min 

f 

J'R-rn 




J{£\c) ct>s[£^,{TZ),TZ]dTZ. 


(17) 


(18) 

(19) 


(20a) 

(20b) 


In the final e xpres sion, T^-min = TZ\c{£miu) and £\c is the value of £ satisfying TZ = 72.ic(£’). We expect the first of the 
two terms in (I20bl) (flux due to diffusion in L) to dominate the loss rate. 
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A quantity that plays an important role in the angular momentum diffusion of orbits near the loss-cone boundary is 


V{E) = 


lim 

K-i-O 


27^ 


lim 


Dnn{E,'R-) 

n 


In this limit, and ignoring diffusion in E, the Fokker-Planck equation becomes 

dN ^ d f dN\ 
d{vt) mz ) 


( 21 ) 


( 22 ) 


showing that D ^ is effectively an orbit-averaged, angular momentum relaxation time at energy E. Another quantity 
that can be expressed in terms of T) is 


qic{E) 


P{E)V{E) 


(23) 


where q\c > 
more useful 


1 defines the full-loss-cone regime. Since the limiting value of 2? as 7?. —>■ 0 may not be well defined it is 
to define V as 


V{E) = 


27^ 


'R-—'R,\c 


DnnjE yTZic) 

TZic 


(24) 


3. DIFFUSION COEFFICIENTS 


3.1. Cohn-Kulsrud diffusion coefficients 

Orbit-averaged diffusion coefficients were derived by iCohn fc KulsrudI (jl978li for stars moving in a Kepler potential. 
Their derivation was based on the theory of random gravitational encounters as developed by Chandrasekhar, Henon, 
Spitzer and others. Cohn & Kulsrud wrote their orbit-averaged diffusion coefficients as 


{AE)t = -Fo + Fi 


{iAEf)t = ^EFo + ^EF4, 


ATZ ATZ 

{AEAn)t = —K4 - —K5, 


(A7^)t = ^ (1 - 27^) Fo + - ~F2 


2E 


E 3 f 


-±F 


2 , ^ Tz 4 7^2 8 7^2 7^2 37^ 

F, + 8-F3 + -—F, + 2 —F, ---Fr 


((A7^)^ = --(l-7^)Fo-2 ^ 


(25) 


where the functions Fi, i = 0 .. .7 are expressed as integrals depending on /. We follow those authors and assume that 
the Fi depend only on the angular-momentum average of /, or 


f{E,t)= [ f{E,n,t)dn. 
Jo 


With this simplification, Cohn & Kulsrud show that the Fi are 

Fo(^)=47rr f~f{£')dS\ 

JO 

F,(£,7^)=47^^^ 'f{E')C,(^,'R^d£', 

where T = AnG'^ml In A, In A is the Coulomb logarithm, and x± = ^ 1 ± (1 — . The functions Ci (s 

have the forms 


^ V (l-x)”/" 


(26) 


(27a) 

(27b) 

E'IE,n) 


with (/, TO, n) integers and Q = {x+ — x){x — X-). Appendix A gives explicit expressions for the Ci and describes how 
they were computed numerically. 

Henceforth the diffusion coefficients (1251) will be referred to as the Cohn-Kulsrud (CK) diffusion coefficients and 
given the subscript “CK.” The subscript t, for orbit-averaging, is understood in everything that follows and will be 
omitted henceforth. 

The quantity ^(f) = [((A7e) )t/{.‘2TZ)]n=nic defined in equation (1241) , which is effectively an orbit-averaged, angular 
momentum relaxation time, can be written in terms of the Fi as 


FiE) = ^ [ 5 Fo (E) + I2F3 (£, 7 ^lc) - 4F7 [E, 7 ^lc)] 


( 28 ) 
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and the quantity q\c defined in equation (E51) is 


qU£) = [5^0 (^) + 12i"3 (£,7^lc) - 4Fr (f,7^Ic)] • 


(29) 


3.2. Resonant relaxation 

The theory of random gravitational encounters that is the basis for the Cohn-Kulsrud diffusion coefficients fails to 
adequately describe the evolution of orbits sufficien tly near to a SBH, where motion is close to Keplerian and where 
orbits maintain their orientations for many periods (jRauch fc Tremainelll996r) . In this regime, it is common to assume 
that random gravitational encounters are still active at changing orbital E and L, according to the diffusion coefficients 
defined above, but that torques due to the nearly-fixed Keplerian orbits are also effective at changing L, sometimes 
on time scales much sh orter than the relaxation time defined in terms of random encounters (jHopman fc Alexander! 

[2?MlEiIorenillM)l . 

We make the same assumptions here, and write the diffusion coefficients that appear in the Fokker-Planck equation 
as (orbit-averaging understood) 


(AS) = (A£’)ck, 

(A7^) = (A7^)cK- 


{{A£f) = {{ASfhK, 


{ASATZ) = (ASATZ) 


(A7^) 


RR 5 


CK, 

{{ATZf) = {{Anf)cK + ((A7^)')RR . 


(30) 


It is understood that the resonant diffusion coefficients describe changes in L in the “incoherent” (as opposed to the 
“coherent”) regime, i.e., on time scales long compared with the coherence time (defined below). 

No very complete theory of incoherent resonant relaxation exists. While the approximate dependence of the diffusion 
rate on energy (i.e. distance from the SBH) is not difficult to derive, until recently, little was known about the angular- 
momentum dependence of the first- and second-order £ -diffusion coefficients. We base what follows on the numerical 
treatment of lHamers. Portegies Zwart &: Merritt! (1201411 . who used an algorithm called TP I (“test-particle integrator”) 
to infer values of the angular momentum diffusion coefficients for test stars orbiting in nuclei with n(r) oc and r~^. 

Those authors expressed the diffusion coefficients in terms of the angular momentum variable 


L 




Lc{E) 

A straightforward transformation yields the relations between diffusion coefficients in £ and in TZ: 

(A7^) = 2£(A£) + {(Aif), {{ATZf) = A£'^{{A£f) 


or 


m= 


1 


2y/TZ 


(A7^) - ^((A7^)2), ((A£)2) = ^((A7^)2) 


Hamers et al. (2014) proposed the following forms for the the first- and second-order diffusion coefficients: 

(A£)rr = Cl A{£) g{£) ((A^)")rr = C 2 A{£) h{£) 


A{a)=al 


Mi,{a) 


M, 


If 


1 Aoh(^) 
N{a) PiaY 

h(^£) = 1 - £ 2 . 


(31) 

(32) 

(33) 

(34a) 

(34b) 

(34c) 


In these expressions, a is related to £ via equation ((S]). The value of was determined numerically to be ^ 1.6, and 
£c was estimated to be ^ 0.7, weakly dependent on a. Hamers et al. suggested Ci ~ C 2 ~ 1. The quantity tcoh(a) 
in equations (l34l) is the “coherence time,” defined as the typical time, for stars of semimajor axis a, to process by an 
angle tt. Hamers et al. assumecQ 

(35) 


+ -1 =+-l _L +-1 

‘•coh — ''coh.M w ''coh.S 


where 


tcoh,M(^) — 


M, 


1 a 


N{a)m- 


-P{a) , tcoh.s(a) = 7 : 7 —-P(a)- 

12 Ta 


(36) 


The latter are averages over eccentricity of the mass- and Schwarzschild apsidal precession times, respectively, assuming 
a “thermal” eccentricity distribution, N{e]a)de oc ede (|Merrittl 120131 §5. 6 .1.1). The Schwarzschilci coherence time is 


^ Some authors, e.g. IHopman fc Alexanderi 112009) and IMadigan et al.l 112011]) . define the coherence time in terms of a dif ference betwee n 
the two precession rates, implying an infinite coherence time at some radius. The reason why this is incorrect is discussed in IMerrittI 1120131) . 
§5.6.1. 
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independent of the mass distribution and is given by 


, , ^ c2a5/2 

tcoh,S\^) — „ s'K/2 

6 {GM,r^ 


, / M, 
5.1 X 10^ , 

VIO 6 M 0 


-3/2 


N 5/2 

- ) yr- 

V 


10“3pc^ 

The mass coherence time depends on the stellar distribution. These two coherence times are equal when 

aN{a) = 12 —-Vg 


(37) 


(38) 


In the approximate theory that motivated the expression (1341) for A(a), the quantity N(a) could mean either “number 
of stars with instantaneous radii less than a” or “number of stars with semimajor axes less than a”. In fact, these two 
functions are quite similar, at least in nuclei with steeply-rising density near the SBH (Appendix B). In his numerical 
experiments, A. Hamers (private communication) was not able to determine which definition of N(a) provided the 
better fit to the data. In what follows, we adopt the first definition, which is simpler to implement in the code: 


N(a) = 4 tt f n{r)r‘^dr (39) 

Jo 

and Mi,[a) = mi,N{a). 

The a— dependence of the diffusion coefficients in equations (|Mll is based on theoretical arguments, but the £— 
dependence is essentially ad hoc. While the diffusion coefficients derived numerically by Hamers et al. were consistent 
with the £ dependence of equations (|34D . the functional forms themselves were not strongly constrained. We now 
consider those functional forms in more detail. 

The variables £ or 72. are bounded, and therefore the flux in 72: 


(j)Tl = -D-Jif — D-jitz 


K 

dn 


(40) 


must go to zero as 72 —>■ 0 or 72 —1. Furthermore this must be true for any /(72) in equation dlH . Hence we require 


Dn — 0 , Dnn 0 

as 72 ^ {0,1}. According to equations (fTSl) and (|32|) . 

Dn = -{An) + ^^((A72)2) = -2i{Ai) + ((Af)^) + £^((A£)2) 
Dnn = \{{An)^) = 2f{iAif). 

Equations (HD) - imply 

(A72)^iA((A72)^), ((A72)^) ^ 0 

i{Ae)^~[e{{Air)], i^{{Aef)^o 

as {72, £} —>■ 0 or 1. 

The expression for ((A£)^) in equation (1341) satisfies these conditions. 

In the case of the first-order coefficient for £, the condition (14^ 1 becomes 

Ci£9ie) = ^-^^m£)], £^{ 0 , 1 } 

i.e. 

£^{0 ,i}. 

Applying this respectively at £ = (0,1} yields 

Co 

Ci = ^ (£ = 0 ), 

Ci( 1 -£^)=C 2 £^ (£ = 1 ) 

and therefore 

£c = ^ ~ 0.577 (C 2 = 2Ci). 


(41) 

(42a) 

(42b) 

(43a) 

(43b) 


(44) 

(45) 

(46a) 

(46b) 

(47) 
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Interestingly, the first-order flux coefficient implied by these functional forms: 

Dn{e) = A{E)' 


-2Ci (^1 - - j + C2 (1 - f) + £C 2 i- 2 i) 
= A{E) [{C2 - 2 Ci)el + { 2 Ci - 3C2Q f] 


(48) 


is identically zero if we choose {C 2 = 2Ci, = 1/3}. Hence, imposing D-ji = 0 at the boundaries implies D-ji = 0 

everywhere - though not necessarily zero flux everywhere. _ 

Based on their numerical experiments, iHamers. Portegies Zwart fc Merritt! (|2014f l suggested a different choice of 
parameters: Ci = C2 and Ic ~ 0.7. Indeed, setting C\ = C2 in equation (I46bl) yields 



£ — — Ki 

'^-V2 

0.707 {C 2 = Cl). 

(49) 

The implied flux coefficient is 

1 

1 



Dn = 

(50) 

zero at .^ = 1 but not at .^ = 0. 





In face of these issues, we returned to the numerical results of Hamers et al. and considered more general functional 
forms for the diffusion coefficients. 

The 77.-diffusion coefficients implied by equations (IMl) have the forms 


(A7^) = A{E) [(2Ci + C 2 ) - (C 2 + 2Ci/£^) 7^] , (51a) 

((A7^)^)=H(^;) X 4C'27^(l-7^), (51b) 


linear in the case of the first-order coefficient and quadratic in the case of the second-order coefficient. Consider the 
more general functional forms 

(A 7 e) = A{E) (a -b -f cE?) , (52a) 

{{/S.nf)=A{E){d + en + fE?) (52b) 

which adds an extra parameter, for a total of six. We now require that these parameters be chosen so as to satisfy all 
the following conditions: 

1. ((A7^)^) > 0 


2. {{ATZf) = 0 for = {0,1} 

3. (A7^) = i4((A7^)2)for7^={0,l} 

It is easy to show that these conditions leave only two independent parameters, allowing the diffusion coefficients (1521) 
to be written as 

(A7^)=H(^;) X C'[l-ry7^^-(ry-2)7^^] , (53a) 

{iAE,f)=A{E) x2C xE,{l-E,). (53b) 

The parameter C is a normalization; based on comparison with equation (1511) . we expect C = 2Ci + C 2 ^ 2. The 
parameter 77 can have any value; it determines the value of TZ at which (ATT.) = 0, i.e. the value of a; = TZq that solves 

{rj — 2) — r]x + I = 0 (54) 

(there is only one root in the range x C [0,1]). r] = 00 corresponds to TTq = 0 and rj = —00 to TZq = 1. 

A. Hamers kindly provided data files with the numerically-determined diffusion coefficients, as presented in the 
Hamers et al. 2014 paper. The best-fitting values of {C, 77 } were determined from these data, at each semimajor axis 
bin, as follows: 

1. Given the numerically-computed diffusion coefficients in £, the 1st- and 2nd-order diffusion coefficients in TZ were 
computed at each of the data points TZi using equations (15^ . To simplify the notation, we call these Di{TZi) 
and D 2 {TZi) respectively. 

2. Data were excluded if they lay outside the range £1 < i < £ 2 , with £2 = 0.95 and £i equal to 1.1 times the 
predicted location of the “Schwarzschild barrier” (as defined below). 

3. The quantity 

^ {Di(7T0 - C [1 - 777 T, + (77 - 2) Uj] }" + [D 2 ( 7 T 0 - 2CTZ, (1 - TT,)]" 

data 

was minimized on a grid in {(7, 77 }. 











0.01 0.1 1 


R 

Fig. 1.— Fits to the diffusion coefficient data from [Hamers. Portcgics Zwart Merri^ l|2Q14l') for the radial bin (a) = 19 .9 m pc. In the 
upper panel, red symbols are — (A7^). Dashed vertical lines delineate the region in which data were fit to the analytic forms ISSll; the latter 
are shown as the solid curves. 


Figured] shows the fits to the data in the bin (a) = 19.9 mpc. 

Figured) shows TZq as a function of a. There is much scatter, but the mean value is close to TZo = 0.5 (77 = 2) and 
there is no obvious trend with a. 

Based on these results, the following functional forms were adopted: 

(A7^)RR = 2A(F;)(l-27^), (55a) 

((A7^)")RR = 4A(F;)7^(l-7^). (55b) 

Note that the TZ^ term in (AT?.) vanishes for the adopted value of 77 = 2. The corresponding diffusion coefficients in £ 
are 

(1 - 

(A£)rr = A(F;)^^^, (56a) 

{iA£f)nK = AiE){l-£^) (56b) 

implying (Ai) = 0 at £ = l/-\/3 « 0.58, consistent with Figure [2 Jd. 

In what follows, equations (|55ll are assumed to define the RR diffusion coefficients in equations ([30|). 

Given these choices, and setting changes in E to zero, the flux coefficients have the simple forms 

Dj^in) = 0, Dnnm = 2 A{E)n (1 - 7?) (57) 


and the 7?— directed flux is 


Fn = -2A[E)n{l-n)^ 


(58) 
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Fig. 2.— The parameter rj in the fits, shown here in terms of the value of TZ where the first-order diffusion coefficient in 'R. is zero (left), 
and the value of i where the first-order diffusion coefficient in £ is zero (right). Dashed horizontal lines are unweighted means of the plotted 
points. 


which has the correct behavior at 7^ = {0,1}. The zero-flux solution is therefore f{TZ) = const, as in the NR case. 
The constant-flux solution, with boundary condition f = 0 at TZ — TZic, is 


f{Tl) = fo log 


1 — TZic TZ 
1 — TZ TZ\c 


(59) 


and the flux is —2A{E)fQ. Since 

f{E)^J= y'/(7^)d7^ = -/olog7^lc, 


(60) 


we can write 

l0g(l/7^1e) 

yielding for the flux 

2A{E)f{E) 
^ "" log(l/7^lc) ■ 


(61) 


(62) 


We expect these expressions to be only approximate since in reality, the steady-state solutions will have nonzero fluxes 
in the E— direction. 

We also give here the expression for V, equation in the RR regime: 


V = 


((A7^)^)t 

27^ 


= 2a1 - 


fM^Y 1 teoh 


- TZic 


\M,J N P2 


(63) 


Using the approximate expression for TZ\c in equation ((TT|) (corresponding to the capture radius for a compact object), 
the quantity gic defined in equation (l23ll becomes 


. X cig / rri-), a 

®c(a) ~ — 1 — ) — ■ 


8 \M, 

Assuming respectively that Uoh = {Uoh.M, Aoh.s} yields 


2 2 / \ 2 
aj m* a ( m* 


N{a) 


Note the interesting result that gic,M does not depend on the mass distribution, while ®c,s does. The former is 


®c,M ~ 0.07 


M. 


^ / M, 


105 


\ 1O5M0 / mpc 


(64) 


(65) 


( 66 ) 















10 


3.3. Anomalous relaxation 

The diffusion coefficients defined in the previous section are affected by general relativity (GR) to the extent that 
GR determines the coherence time; the latter defined as a typical precession time of all stars at a given radius. 
Another GR-related phenomenon is the “Schwarzschild barrier” (SB), the tendency of single, very eccen tric orbits not 
to di ffuse in L below some definite value at each a. The SB was first observed in A^-body simulations (|Merritt et al.l 
[Ml , as a locus in the logo vs. log(l — e) plane where the orbital trajectories “bounced” in the course of their 
RR-driven random walks in L. The same study revealed that orbits experiencing the “bounce” were of such high 
eccentri city that their GR precession times were short compared with those of typical (i.e., higher-L) stars at the 
same a. iHamers. Porteeies Zwart fc Merritd (j2014[ l coined the term “anomalous relaxation” to describe the behavior 
of orbits in this high-eccentricity regime. 

Two analytic expressions have been proposed for the location of the SB. The first is based on a comparison of the 
GR precession time with the time for the VN torques to change L (|Merritt et al.ll2Ml : 






M. 


aJ 


Nia). 


(67) 


The second (jHamers, Portegies Zwart fc MerrTttll2014ll compares the GR precession time with the coherence time: 




^coh(^^) 

' 1 ^ 


( 68 ) 


In spite of their disparate functional forms, the two expressions yield numerically similar relations for TZsb{o-) in a 
wide range of nuclear models . However the former relation appears to more accurately rep roduce the barrier location 
in numerical studies to date (iMerritt et al.ll2Ml [Hamers. Portegies Zwart fc Merritt! [2M1 . 

Appendix G deriv es expressions for the diffusion coefficients in this regime, based on a simple Hamiltonian model 
(jMerritt et al.ll2??Tl : 

2£3 _ ei 


{M) ^ —, {{MY) ^ - 


(69) 


where r(a) = tcoh(a)/(A^)^ and 


A 


Vn — 


Mi, {a) 


2./N{Y M, Tg 


(70) 


[Hamers, Portegies Zwart fc MerrittI (|2014ll verified the ^dependence predicted by equation (IMl) via numerical experi¬ 
ments. Expressed in terms of TZ via equations (1321) . the diffusion coefficients become 


(A7^)AR = ((A7^)")AR = 

T T 


(71) 


Equations (ED replace equations (1551) when TZ < TZsb- We note that the values of the numerical coefficients in these 
expressions have not been well determined by the numerical experiments to date, and the simple model in Appendix 
C is not likely to be valid when TZ is so small that the GR precession time is much shorter than the coherence time. 
Additional A-body experiments are needed to elucidate the functional forms of the anomalous diffusion coefficients. 
Soluti ons to the Fokk er-Planck equation containing these diffusion coefficients will be postponed to Paper HI in this 
series ()Merrittil2015bll . 


3.4. Gravitational radiation 


First-order diffusion coefficients in £ and TZ also exist that describe changes due to emission of gravitational waves. 
At the lowest (2.5 PN) order, the orbit-averaged rates of change of a and e are: 


(Aa) = 
(Ae) = 


64 G^MYrOi, , 73 2 , 37 4 ^ 

Tc5a3(l-e2)"/2 j ’ 

304 G^MYrUi^e / 121 Y 

13 (1 — e^Y^"^ k 304 J 


(72a) 

(72b) 


(equations 4.234 from IMerritt! (j2013l) after replacing mim in those expressions by MY). Expressed in terms of 
£ = GM,/{2a) and 7^ = 1 — e^. 


(Af) 

(A7^) 


2£’2 

GM, 


(Aa) = 


-2\/r77^(Ae) 


2720 1 ^ ^ ^ X 

3 M, GM,c^ 7^7/2 V 425 425 J ’ 

2720 1 / 540 42 I 

3 M, GM,c^ 7^3/2 425 425 J ' 


(73a) 

(73b) 
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4. NUMERICAL ALGORITHM 
4.1. Choice of grid and units 

Solutions are obtained on a {Nx x Nz) grid in (X, Z), where 

T n 2 


X = In ii = In 


= In (l - e^) , 


[Lc{£). 

Z = \n{l +I3£*) =ln{l +IdS/c"^) 

(Figure [3]). Typically Nx = Nz = 64. Here and below, starred variables are dimensionless, based on code units which 


(74a) 

(74b) 


are 


G = M. = c=l. 


(75) 


Thus £ = c^£* and L = {GM,/c)L*. The angular momentum of a circular orbit, Lc = \/GM,/{2£), becomes in 
dimensionless units L* = l/{2£*). The squared angular momentum of an orbit with periapsis at r = 8rg is 


RiJe*) = 32£* (1 - 8 £*), £* < — . 

16 


(76) 


In the code, all diffusion coefficients are divided by the factor T = In A. The dimensional factor /q multi¬ 

plying /* is taken to be 


/o — 


These two choices determine the unit of time as (T/o) ^ or 


9 

-1 


M, 

TO* 


1 


[t] =to = 


M, 

47 r In A \ to* 


2-7 


AG \ ^ \ ('fUl.'] - 1 

TO* / \ r„ J \ c J 


= 1.7 X 10"V 


/InA 


-1 


M,/m 

105 


3-7 


(77) 


(78) 


In these expressions, and 7 are free parameters, but they have a simple physical interpretation in the case of a 
cluster with a power-law density: 


i(r) = no ( — 


no = 


3 — 7 M, 1 
2tt to* 


(79) 


for which is the radius containing a mass in stars equal to 2M,. The isotropic / corresponding to the density GHl) 
is 


f{£) = for{£*) = /oC'( 7 )f*"-'/^ G( 7 ) = ■ 


(80) 


Initial conditions were typically chosen to be equations (1731) - (1801) or some modification, e.g., / might be set initially 
to zero inside the loss cone, as described in more detail below. When the initial / was constructed in this way, Cm and 
7 retain their meaning as the parameters describing the unmodified power-law model. 

The initial density normalization (stars per unit volume) is fixed by the choice of = rzuffg. At any later time, 
the dimensionless number density: 

P* ^3^,* 


is related to the true number density by 


n*{r*) = J rd^v* 

(r) = [ fd\ = cVo f rd\* = 1 ^ 

J J fqX'^gJ 


Similarly the number of stars within r is 


X(< r) = 


7-3 


V™* 


N*{< r*) 


(81) 


(82) 


(83) 


where N*(r*) = 47 r n*(r*)r*^dr*. These expressions are used when evaluating the diffusion coefficients in the RR 

and AR regimes. Combining equations (1781) and (1531) . the rate of loss of stars from within a sphere of radius r is given 
in terms of dimensionless quantities as 


=47rlnA- 


dt 


dt* 


(84) 
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CO 
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+ 



c 


M 



-15 -10 -5 0 


X = In i? 


(b) . 

• 

•1 

• 

• 

.bi 

—^ • 

• 

(c) 

I' 

• 

1 

• 

.^i+i,j+i 

0 1 ,j 

\ .bi 


o^i “1 

1 okj-F 



Fig. 3. — (a) Solution grid; dots denote cell centers where the fij values are defined. Fluxes are defined on the cell boundaries. Thick 
solid curve is the boundary of the loss cone, 7^ = TZidS). (b) Schematic diagram showing the fluxes associated with grid cell (c) 

Finite-differencing for grid cells adjacent to (above) the loss cone boundary is affected by the lack of knowledge of / inside the loss cone. In 
the case of “empty loss cone” boundary conditions, shown here, / is assumed to be zero inside the loss cone (open circles). Cohn-Kulsrud 
boundary conditions are expressed purely in terms of the values of / near the loss-cone boundary; no derivatives are evaluated there and 
the values of / inside the loss cone are not needed. 


The dimensionless function gic(^), equation (|23p . that appears in the Cohn-Kulsrud boundary layer treatment also 
depends on Since 

Vi£)=rfoV*{£*), (85) 


we can write 


where 


qic{£*) = 2V2 tt‘^ 


7-3 


V*{£*) 1 
£:*3/2 ^ 


V*{£*) = 


^nn ^ic) 

'R-ic 


( 86 ) 


(87) 


is the dimensionless form of the function defined in equation (1241) . For instance, in the case that diffusion in L is due 
only to non-resonant relaxation, equations (IMl) and (|29ll imply 




2-\/27r^ TO* 

3 Jf, 


InA — 


7-3 


5Fo* -f I 2 F 3 * - 4Ff 
£*^/'^n\c 


( 88 ) 


4.2. Basic numerical approach 

As in iCohn &: KulsrudI (I1978I1 . the distribution function / is evaluated at the cell centers while the flux coefficients 
and fluxes are evaluated at the cell boundaries (Figure [3li). The rate of change of J'/ in a cell is then equated to the 
net flux through the cell boundaries (Figure [3})). The dimensionless evolution equation: 


dt* 




(89) 
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was written in implicit discrete form, with the time derivative represented as a backwards difference and the flux 
coefficients evaluated at the previous time step. In discrete form, the Nx x Nz difference equations become 


/, 


(n+l) 


_ 

J a 


At 


i+l 

k—i — 1 l—j — 1 


i+l 


(90) 


where = f{Zi,Xj) evaluated at the nth time step. Solution of the matrix equation (TOl) is carried out using 
the fortran program f07acf from Numerical Algorithms Group. This routine uses iterative refinement to produce a 
solution with double-precision accuracy. Iteration was found to be unstable when the time step, At, was taken to be 
too large; the maximum allowable At was found to decrease with increasing number of grid cells. Less accurate matrix 
i nversion routines were fo und to be faster but also more susceptible to instability. 

iCohn fc KulsrudI (119781 1 employed a linear grid in TZ extending from 7?. = Oto77. = l. As a result, the dependence 
of / on 77. near the loss cone was not well resolved. The use here of a logarithmic grid in 77 yields many more grid 
points at small 77, permitting a more careful treatment of the solution near the loss cone, as discussed below. 


4.3. Boundary conditions 

Two choices for the the boundary condition at small 77 are implemented. The simplest is an “empty loss cone” 
(ELC) condition: f{£,TZ) = 0 for 77 < TZic{£). To implement this condition, cells that contain the curve 77 = 77ic(£’) 
and for which the cell center is outside the loss cone (“loss-cone cells”) are identified, and for these cells, a record is 
made of the adjacent cells for which the cell center is inside the loss cone, i.e. for which TZij< TZ{£ij). There are 
nine possible combinations of adjacent cells for which this condition can be satisfied; Figure illustrates the most 
common, consisting of three adjacent cells. The finite-difference expressions for the flux derivatives in the loss-cone 
cells are written with the appropriate /ij-values set to zero. 

The second choice of small-77 boundary condition was based on the Cohn-Kulsrud (1978) boundary layer solution. 
In the spirit of that derivation, the evolution equations near the loss cone boundary are first simplified by ignoring the 
contribution of gradients in £ to the flux, which allows equations (II3 to be written as 


df df 

4>e = —Dsni^ — Dsf, (j)Tz = —D-jm-— — D-jif. 

oTZ oK 

In the Cohn-Kulsrud solution, the 77-directed flux per unit of £ across the loss cone boundary: 

Fi£) = -£ dU = J{£) [(j)n{n = 1) - <(.7^(77 = 77ie)] = -77(f)07^(7^ic) 


(91) 


(92) 


is given by 


F{£) = iTT^Ll{£) 77le(f) /(f,77le) a<lic) 


(93) 


where 


i{x) = 1 - 4 ^ 




(94) 


and the am are the consecutive zeros of the Bessel function Jo{a) (|Merritt]l2013L equations 6.58-6.62). Thus 


= -^-—tzU£) fU£) 

TT LrlVlt 

where fic{£) = f{£,TZic). The gradient in / at 77 = 77ic in the Cohn-Kulsrud solution is 

V 977 /i^ qicTlic 

which is inserted into the first of equations (ITO to give the flux in £: 

^e{£,Tl,,) = -Dsn{£,ni,)^^^ - De{£,'Jli,)f,,{£). 

qic 77ic(t) 


(95) 


(96) 


(97) 


Equations ((55|) and (IWl) are the adopted boundary conditions. These expressions are applied at the inner edges of any 
loss-cone cell, while the general expressions m for the flux are applied on the outer edges of those cells. 

The summation in equation (1941) is slowly converging when q is small (Figure 0^). An approximate expression is 


= 


(cc^ -I- 


(98) 


for which the relative error is ~ 1% at all g < 1 and tends to zero for large gic (Figure|4|D). Equation (IMl) was adopted 
for the code. 
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Fig. 4.— (a) Convergence of the series in equation Jolt . Each curve is labelled by rumax, the number of terms included in the summation; 
the ordinate is the relative error, (b) So lid curve shows ^{q) computed from equation (I94II using 10^ terms. Dashed curve is the absolute 
error in the approximate expression ll98l l. 

Note that the values corresponding to grid cells inside the loss cone are not used and therefore are not advanced 
in time. Indeed in this region, / is not expected to satisfy Jeans’s theorem and a proper treatment would require 
inclusion of r— as an additional independent variable. 

The boundary condition at 7?. = 1 is (jin = 0- 

Two possibilities were implemented for the boundary condit ion at small 8, i.e. far from the SBH. The first consists 
of fixing /(fmin, 7?.) at its initial value. ICohn fc KulsrudI (|1978n implemented a similar boundary condition. This choice 
implies a nonzero flux of stars into the region of integration, although that flux can be very small if fmin is small. The 
second choice consists of setting to zero the £— directed flux at £ = fmin, i-e. requiring 


df df 


(99) 


This boundary condition is more in keeping with iV-body simulations that have a fixed number of stars. The “zero-flux” 
boundary condition is implemented as follows. (1) At each time step, inversion of the matrix Aijki is first carried out 
specifying no change at f = fmin, be., Z = Z^z- (2) Equation (IMl) is then used to solve for /{Znz, Xj),j = 1,..., NX. 
In this step, derivatives are expressed in terms of the unknown values of / at Zpfz and the just-computed values at 
Zffz-h allowing solution of / by inversion of a tri-diagonal matrix. 

4.4. Initial conditions 

The initial f{£,TZ) was typically based on an isotropic power-law model, equation (1501) . but with a modified 72.- 
dependence to account for the presence of the loss cone. The simplest choice is to set / = 0 for 72 < 72ic. A more 
natural choice is to set the initial / to 


f{£,n) = /{£,!) 

= 0 


ln(72/72ic) 


ln(l/72ie) ’ 
72 < 72ic 


72 > 72ic 


( 100 ) 


which is the approximate (small-72) steady-state solution for an empty loss cone. Other choices for the initial conditions 
are described below. 
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5. EXAMPLES 

Subsequent papers in this series will describe the evolution of f(E,L) due to the various physical mechanisms 
represented above. Here we use the algorithm to explore some simple dependencies of the classical Bahcall-Wolf 
solution on and on time. 


5.1. Models with classical dijfusion coefficients; empty loss cone 

The number of parameters that define an integration can be minimized by assuming empty-loss-cone (ELC) boundary 
conditions, and adopting the classical (CK) diffusion coefficients. The stellar mass, m*, then appears only in expressions 
like (1751) that set the scaling of the time or density, and any initial f{£,TZ) should evolve to the same equilibrium / 
modulo a known scaling in amplitude. Furthermore this steady state should be similar to the Bahcall-Wolf “zero-flux” 
model, n ^ f ^ the differences being due to the fact that those authors assumed / = f{£,t) and allowed 

loss of stars to the SBH only via diffusion in energy. 

The code was tested under these assumptions, starting from power-law models like those of equation (1801) and various 
values of 7 . The isotropic models were modified initially as in equation (llOOl) so that / fell to zero at 72. = 72ic; the 
ELC boundary conditions then guaranteed that /(72ic) remained zero at all 72 < 72ic. /(£min,72) was fixed to its initial 
form. 

It is permissible and convenient to scale these models by assuming that the final mass density has a specified value 
at a specified radius. Let that radius be arg and let the density at that radius be p(arg). It may be shown, using 
equation (fTSl) , that the elapsed time is then given in terms of the dimensionless time t* by 


t= ■ 


1 


47 r In A 
!7.4 X 10^^ 


M, 

m* 


(?) 


/InA 


-1 


P{(^rg)rl 

M, 


n*{a) t* 


f M, 


-1 


105 J MqJ [lO6M0pc-3_ 


p{arg) 


n*[a) t* yr 


( 101 ) 


where n*(a) is the final, dimensionless model density (cf. equation IHTI) at r = avg, i.e. at r* = a. The loss rate, 
equation (ISil) . can similarly be related to the dimensionless loss rate as 


—— = 47 rln A 
dt 



p{arg)r^ ^ 1 dN* 

n*(a )2 dt* 


a; 1.49 X 10"^® 



M, 

IO 6 M 0 


Pjo^rg) 

10^Mqpc~^ 


1 dN* 
n*(a)^ dt* 


yr 


( 102 ) 




Fig. 5.— Evolution due to classical relaxation of two models assuming empty-loss-cone (ELC) boundary conditions. The initial model had 
n ~ r~^ (left) and n ~ (right), shown as the dashed lines. Other curves are at times (1, 2, 3, 4, 5) x 10^ yr (left) and (0.01, 0.1, 0.5, 7) x 

10^ yr (right) based on the scaling described in the text; line thickness increases with time. Dotted lines show dlogp/dlogr = —7/4, the 
Bahcall-Wolf slope. 


Figure [5] shows the evolution of p(r) for integrations with 7=1 and 7 = 9/4. The units of time and density were 
fixed by assuming a final density at lO^r^ of IO^Mqpc”^; the other parameters were set to the fiducial values in 
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equations (IIOII) and (I102p . e.g. M, = 10® Mq, hence 10®rg ~ 4.8 x 10 ^ pc; the radius containing a mass in stars of 
2M, at the final time step, the “gravitational influence radius,” is ~ 3 x 10®rg, roughly equal to the outer radius of 
the solution grid. In the model with initially shallower slope (7 = 1 ), the density evolves “from the outside in” toward 
the steady state, while in the model with 7 = 9/4 the density evolves at all radii at roughly the same rate. 

The final density profile is the same in these two models, as expected, and is close to a power law: 

p{r) - r-\ 1.65 <5< 1.70 (103) 

at lO'^rg < r < lO^rg. This is slightly shallower than the canonical p ^ of the Bahcall-Wolf solution. The 

difference can be attributed to the depletion of orbits in the loss cone, the effects of which are progressively more 
severe at energies close to the SBH. 




time (yr) 


time (yr) 


Fig. 6. — Evolution of the peak density (left) and of the total loss rate (right) for three integrations assuming classical diffusion coefficients 
and ELC boundary conditions. The initial density slope was 7 = 1 (triangles), 7 = 1.5 (open circles), and 7 = 2 (filled circles). Scaling 
assumes a final density at 10® of 10® Mq pc“®. 


Figure m shows the evolution of the peak central) density, as well as the total loss rate TV, for three integrations 
starting from 7 = {1, 3/2, 2}. In each case a steady state is reached in a time of ^ 5 x 10® yr. This can be compared 
to the classical relaxation time, a standard expression for which is 

Tr{r) = , (104) 


G^mp In A 


0.95 X 10^° 


P 




yr 


. 200 kms-V VlO®M 0 pc-V \Mq) \ lb ) 

(jMerrittll2013L Eq. 3.1). Evaluating equation (11041) in the final model, assuming n{r) ^ m* = 10“®M, = IOMq, 

In A = 15, and 


cr(r) 


1 GM, 


1 + 7 r 


1/2 


183 


10 ®r„ 


- 1/2 


km s 


-1 


yields 


Tr{r) « 7 X 10® 

consistent with the observed equilibration time. 


10 ®r„ 


0.20 


yr, 


(105) 

(106) 


5.2. Models with classical diffusion coefficients; Cohn-Kulsrud loss cone 

The simple rescaling defined in the previous section no longer works when adopting the Cohn-Kulsrud (CK) loss-cone 
boundary conditions, equations (IMl) - (1^71) . since the dimensionless function qyfiSfi) depends separately on the factor 

In A 

M, 


(107) 
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(equations uni US]) . Two models with different initial f’s will not necessarily arrive at steady states that are rescaled 
versions of one another unless the final qidS) happens also to be the same. 




I°gi0 '■/'■g 




Fig. 7.— Final quantities for a set of integrations adopting the Cohn-Kulsrud (CK) boundary conditions, classical diffusion coefficients, 
and various values of Thick solid (black): qic = 0 (empty loss cone; (b)-(d) only); dashed (blue): rrii^lMm = 10“^; dot-dashed (red): 

10“^; dotted (purple): 10“^; dash-dot-dot (orange): 10“"^; thin solid (black): 10“^. Open circles mark the energy where qic = jlnT^icI? 
roughly the expected energy of transition between the empty- and full-loss-cone regimes. 


Differences in scaling will be minimized if the initial model is close to the final model. In this section, initial conditions 
are chosen to be close to the Bahcall-Wolf form, / ~ Differences in the final state will be due to differences in 
the value of (I107L i.e. to differing degrees of loss-cone “fullness”. Parameters common to all the integrations in this 
section were 

7 = 7/4, — = 10^ lnA=15. (108) 

Various values were adopted for of rrii^/M,, implying different forms for qic{£,t = 0). The initial density profile is 
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approximately 


p{r)'. 


3 — 7 M, f r 


2tt \rm 


!2.0 X 10- 


and the unit of time becomes 


/ M, 

\10^MqJ VlOpc. 

' M,/m^ 


-3 / \ -7/4 

r ' 


[t] ft! 3.0 X 10^ 


105 


— ) 

I 71 


yr. 


These initial models are expected to reach steady states characterized by 

7^ 


for qic 1 (“empty loss cone”) and 


/(I) 1 

^ ln(l/7^lc) V^ic 
/(72.) ~ const. 


(109) 

( 110 ) 


( 111 ) 

( 112 ) 


for Qic » 1 (“full loss cone”). 

Figure [3 plots four dimensionless quantities associated with the steady-state solutions. 

(a) qic{£)- The transition from empty- (qic 1) to full- {qic ^ 1) loss-cone regimes takes place at progressively larger 
binding energies, i.e. smaller radii, as m* is increased. For the smallest value of m* adopted here (m*/M, = 10“^) 
the final model is essentially in the empty-loss-cone regime everywhere. 

(b) The logarithmic derivative of the configuration-space density, dlogn/dlogr. At intermediate radii, 10^ < r/rg < 
10^, the slope is close to —1.7 in all final models, slightly shallower than in the (isotropic, scale-free) Bahcall-Wolf 
solntion. This is due to the progressive depletion of the loss cone near the SBH, as in the ETC examples above. That 
depletion becomes less severe as m* and hence qic are increased and the slope in the models with the largest m* 
approaches most closely to —7/4. 

(c) The magnitude of the 77.-directed flux, (/-r, = -D-jinidf/dlV) — D-pif, evaluated at 77. = TZic{£). The empty-loss- 

cone model exhibits an approximate power-law dependence, (j>Ti,ic ^ at small £* (large radius). In the models 

with finite rrii,, the flux drops sharply beyond the energy where qidS) ln77ic. 

(d) The quantity J{£)4>tz,\c{£)^ whose integral d£ is proportional to the integrated loss rate. Since J oc £“5/2^ 
•J^ £~^ for small £ in the empty-loss-cone model, implying an integrated loss rate that diverges as £~^ ^ r. 
When m* is non-zero, £f(j)-R,\c instead “levels ont” roughly where gic ^ \ ln77ic|, implying a finite total loss rate. 

Equation (I5H) . together with the parameters (11081) . yields a relation between the dimensional and dimensionless loss 
rates: 


dN{< r) 
dt 


= 47rln A- 


dN*{< r*) 

\rg J dt* 


3.8 X 10 


-14 (_ 


M, 


V lO^M© J dt 


dN* 


yr 


(113) 


Figure|5^ plots this function, at the final time, for the five models with non-zero gic, assuming M, = lO^M©. We can 
comp are these loss rates to those predicted by a simple model based on the assumption of an empty loss cone (jMerrittl 
I2013L equation 6.91): 


N{< r) 


Jn. 

InTT-ici Tr{r) 


3 (l-b7)3/2(3-7)^ jV(<r^)2Gi/^w;inA 
477 (9/2-27) M.3/V^/^ln77ic 

3^ InA r 

ln77ic y 


5.6 X lO”'^ 


M, 

IO^Mq 




(114a) 

(114b) 

(114c) 

(114d) 


(The second line follows from equations (110411 and (110911 : the third line sets 7 = 7/4; and the final line uses the 
parameters (11081) and appr oximates 77ic by ridrm = Recalling that rm/rg = 10® for the models of Figure|Sl 

we see that equation (11141) correctly predicts the loss rates in the Fokker-Planck models for the case of small to*/M,, 
i.e. in the q\c —?> 0 limit. In the case of an empty loss cone, there is a (linear) divergence of the integrated loss rate 
with r; when m* is finite, there is a transition to the q\c > 1 regime in which the loss rate is given approximately by 


N {< r ) KiA'K ( 77ic(r) dr 

dreHt Py^) 


( 115 ) 
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Fig. 8.— Left panel. Dimensional loss rate, defined as number of stars per year scattered into the SBH from radii smaller than 
r, for the five models of Figure [3 “Radius” here means GM*/(2£^) and M* = 10 ®Mq has been assumed. Line colors and types 
have the same meanings as in Figure [7| Circles indicate the radii where q\c = |ln7^ic|. Right panel. 7^-dependence of / in the 
steady-state model from Figure [7| with = 10“® (NR diffusion coefficients). Circles are f('R.;S) at five values of S, for which 

Q'ic(^) ~ {45.4,5.44,0.816,6.89 x 10~®,6.14 x 10~®}. Curves are equations JllTI t. after normalizing such that the numbe r of stars at each 
S is the same in the numerical and analytic models. The vertical tick-marks at three values of 'R. show 7^o(^)? equation (I118I I. 

where Tcrit is the smallest radius for which qic exceeds one (IMerrittll2013L equation 6.91). For these models, equation 
(I115|l predicts a contribution from the full-loss-cone regimes that scales with radius as 

N{< r) oc r > rent (116) 

consistent with the “leveling out” of the curves in Figure H] at radii greater than ^ rcrit- 
It is interesting to com pare the 7^-dependence of / i n the final models with the form that is often assumed when 
computing loss rates le.g. iMagorrian fc Tremaind 1 199811 : 


/(7^) = f{TZie) + In (7^/7^lc), Teic < < 1, 


= ■ 


ln(l/7^lc) 

/(I) 


l + 9ic ^(9)ln(l/7^lc) 


(117a) 

(117b) 


Equation (I117a|) is a steady-state (constant f lux) s olution of the diffusion equation in TZ if the low-72. form s of the 
diffus ion coefficients are used, and equation (Ill7bll follows from the Cohn-Kulsrud treatment (as given in iMerrittl 
I2013L equations (6.59), (6.61)). Figure |BI) plots /(72), at various £, in the final model from the integration with 
rrii^/M, = 10“®. Equations (11171) are overplotted, after normalizing to give the same total number at each £. The 
agreement is reasonably good, verifying that the boundary conditions have been correctly implemented. But the 
dependence of / on 72 near 72 = 1 deviates from the simple logarithmic form of equation (I117al) , due to the fact that 
the correct - not small-72 asymptotic - forms of the diffusion coefficients are used in the numerical code. 

In the case of a non-empty loss cone, the / = 0 intercept of the “external” (72 > 72ic) solution occurs at 72 = 72o < 72ic, 
where 

72o(<?ic) =72ie(£:)e-«''/«(«'=) (118) 

(|Merrittll2013L equation 6.65). In FigurelHla, 72o as given by equation (IllSp is indicated by vertical marks, and appears 
quite consistent with the value of 72 at which the numerical solutions are tending to zero. 


5.3. Models with classical and resonant diffusion coefficients 

The form of steady-state solutions derived in the preceding two sections depended on the value assumed for the 
stellar mass m*, insofar as m* determines qic- But the dependence was found to be weak, and all of the steady-state 
solutions could be rescaled in / and £ to approximately the same f{£). 

This scale-free property is lost if the resonant relaxation diffusion terms are included. Sufficiently near the SBH, 
changes in angular momentum will be dominated by resonant relaxation while changes in energy will still be dominated 
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by the classical diffusion coefficients. A characteristic radius appears: the distance, req, from the SBH at which the 
timescale for changes in L due to resonant relaxation is the same as the timescale for changes in L due to classical 
relaxation. That radius is given roughly by the solution to 


mi,N{< Teq) = AnrM,, 


Anr 


0.68 1 
(3 - 7)(1 + 7 ) 3/2 log^ 


(119) 


(|Merrittll2013L Equation 5.243). Here N{< r) is the number of stars within a sphere of radius r and t^oh = ^coh.M 
has been assumed; the expression for Anr additionally assumes n{r) oc r~'^. Setting log A = 15, the enclosed (stellar) 
mass at req works out to be ^ 10“^M, for all 0.5 < 7 < 2.0. 

Consider a star whose energy diffuses below Egq ^ —GM,/req- The timescale for diffusion in L will suddenly 
decrease, and the star will be lost to the SBH, in a time much less than the characteristic time for changes in E. A 
steady state will eventually be reached, but only after \df /dE\ at E ^ E^q becomes large enough to drive a flux in E 
(due to classical relaxation) that equals the flux in L (due to resonant relaxation). The expected result is a depletion 
in / at A ^ Aeq and a low configuration-space density at r < req - a “core”. 




Fig. 9.— Evolution of p(r) in two models with the same initial density, p(r) oc r~^. Left: classical diffusion coefficients were used. Right: 
resonant relaxation terms were included as well. Other features of the integrations are specified in the text; the adopted parameters are 
appropriate for the nuclear cluster of the Milky Way, with M* = 4.0 x 10 ®Mq. Densities are plotted at times t = {0, 2,4, 6, 8,10} x 10^ yr. 
The t = 0 model is shown by dashed lines, and line width increases with time. Dotted lines show the Bahcall-Wolf slope, dlogp/dlogr = 
—7/4, and the slope corresponding to a sharply truncated f{S), dlogp/dlogr = —1/2. Both models have an integrated stellar mass, 
within one parsec, of approximately 1.0 x 10®at t = 5 x 10^ yr, consistent with the mass measured in the Milky Way nuclear cluster 
aSchodel et al.lf^^ . 


Figure [9] shows the evolution of p{r) in two integrations: the first using the classical diffusion coefficients, the second 
including the resonant-relaxation coefficients as well. Both integrations adopted parameters appropriate for the nuclear 
cluster of the Milky Way. The quantity (m*/M,) In A that appears in equation (1551) for qi^: 


rrii, 

W, 


InA = 3.75 X 10“® 



M, 


4 X IO 6 M 0 



( 120 ) 


was set to 3.75 x 10 


and the radius of the loss sphere, which for a main-sequence star is the tidal disruption radius: 


He 


. R* 
0.7AU 

Rq 


I'M, 1 

4 X 10®/ 


( 121 ) 


was set to 0.7 AU. These choices were motivated by the fact that the main-sequence turnoff mass in the Galactic center 
is ~ IM 0 , and the red giants that a re believed to dominate the number counts of the “late-type” stars probably have 
roughly this mass (|Dale et al.ll2009fl . (Due to the relative shortness of the red giant evolutionary phase, most stars are 
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expected to be disrupted while still on or near the main sequence (iMacLeod et al.ir2012fl l. These choices determined 
the location and flux at the loss-cone boundary via the expressions (EH), (EZl). 

The initial density normalization was chosen, after some trial and error, to give an integrated mass within 1 pc of 
^ 1.0 X at t = 5 X 10® yr. This is roughly the value inferred from dynamical analyses of stellar velocities 

(|Sch6del et al.11200^ . 

Figure [9] shows that the integration based on the classical diffusion coefficients nearly reaches the p oc r Bahcall- 
Wolf form after 10 Gyr. Deviations from that form are apparent inside ~ 0.01 pc even at this late time; however seen 
from the Earth, that radius would subtend an angle of only ^ 0."25. 

Inclusion of the resonant-relaxation diffusion coefficients implies a different steady state. Inside ^ 0.05 pc, the density 
profile remains much shallower than the Bahcall-Wolf form, with p ^ r~^'^ - the functional form that corresponds to 
an / that is fully depleted at high binding energies. The radius of this “core” is consistent with the prediction made 
above O.OIcm) given that Ri 2.5 pc. 

Since the radius of the core is a function of the density normalization, steady state solutions in the presence of 
resonant relaxation are expected to depend on (at least) one more parameter than i n the classical case. A more 
general exploration of solutions like these will be presented in Paper 11 from this series (|Merrittll2015aH . 


I thank A. Hamers for providing data from his TPI code that was used in deriving the functional forms of the RR 
diffusion coefficients in 113.21 and H. Cohn for helpful discussions. This work was supported by the National Science 
Foundation under grant no. AST 1211602 and by the National Aeronautics and Space Administration under grant 
no. NNX13AG92G. 

APPENDIX 

Numerical Evaluation of the C, 


The functions Ci{s = £'/£,TZ) that appear in equation (I27bl) . 

pS fx- 


F, {£,n)= 47rr 

are given by iCohn fc KulsrudI (|1978ll as 

Ci = — J dx Q 

C 2 = — [ dx Q 


£' 


f{£')a{j,n]d£' 


- 1/2 X (1 - 

(l-Cc)l/2 ’ 

- 1/2 

(l-x)3/2 ’ 

C3 = - /■ dx Q-1/® 

' (1 - x)l/2 ’ 

- 1/2 (1 - Sxf/^ 
(l-x)l/2 ’ 

- 1/2 X (1 - Sx ) 3/2 
(l-x)3/2 ’ 

- 1/2 - sxf^'^ 

(l-x)5/2 ’ 
-1/2 ^3(1 - Sx)3/2 
(l-x)3/2 ■ 


Ci = - 


dx Q 


Cs = — / dx Q 


Cr = - 


dx Q 


C 7 = — / dx Q 


(Al) 


In these expressions, Q = (x+ — x)(x — X-) with x± = (1/2) [l ± , and the limits of integration are X- and 

max [x-, min (x+, s“^)]. These functions are independent of / and so can be evaluated on a fixed numerical grid. The 
grid axes were chosen to be {72., W} where 


W = 


2T-1 + 




and T = s ^ = £!£'. Thus 0<1T<1, 0<72<1, and 


T = 


i 1 - \/r^^-1-VF (i-p 


(A2) 


(A3) 


Typically the number of grid points was 256 x 256. 
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The Ci are finite for all {TZ, W} with the exception of C 2 and Cg, which diverge at 

7^ = 0, W = T = s=l. 


Along the TZ = 0 border, both functions diverge as (1 — W) as IT —>■ 1, while along the W 
as TZ~^I'^ as 7^ —>■ 0. To deal with the divergence, C 2 and Cg were multiplied by the function 


V} + 256 (1 - W f 



(A4) 


1 border they diverge 


(A5) 


before storing their computed values on the grid. Integrations were carried out with routine dOlapf from the Numerical 
Algorithms Group fortran subroutine library. Accuracy of the integrations was checked by comparison with analytic 
expressions that obtain along the grid boundaries; for instance, for TZ = 1, Ci = gi{s) where 

51 =52 = V2 - s, gs = - s, 

54 = 55 = 56 = (2 - sf^^ , 57 = i (2 - . (A6) 

Once the values of the Ci had computed on the {TZ, W} grid, the NAG routine eOldaf was used to fit a bicubic 
interpolating spline to the computed values. During integrations of the Fokker-Planck equation, values of the Ci 
between the {TZ, W} grid points were then computed from the spline coefficients. 


APPENDIX 
Nr (a) vs. Na{a) 

This appendix compares two quantities: 

1. Nr{< a), the number of stars with instantaneous radii less than a; 

2. Na{< a), the number of stars with semimajor axes less than a. 

A power-law dependence of number density on distance from the SBH is assumed: 


n(r) = no 


ro 


ipir) = 


CM, 


so that the number of stars instantaneously below r is 


Nr{< r) = 


Att 

o - noro I — 

3-7 V’"o 


3-7 


The distribution function is assumed to be isotropic; Eddington’s formula gives 


fi£) = fo 


5 = 7-3/2 


1 r(7-H) nprg 07-3/2 

(274)3/2 r(7-1/2) 

The number of stars per unit of binding energy is 


N{S)d£ = ATr‘^p{£)f{£) d£, pi£) = 


so that 


N{£)d£ = 


773/2 


2 r(7-l/2) 

and the number of stars with binding energies greater than £ is 


iV,(> £) = n,rS 


Setting £ = CM,/{2a) yields the number of stars with semimajor axes less than a: 

22-7773/2 r (7 + i) 

Na{< a) = 


(3-7) r(7-1/2)^^°^° V^'o 


(Bl) 

(B2) 

(B3a) 

(B3b) 


p{£) = (gm.)3£:-5/2 

(B4) 

{GM,f~^ norl £^-^d£ 

(B5) 


(B6) 


(B7) 
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Equations (IB2I) and (lB7l) can be written 


Nr{< a) = Ni{^) norp 


ro 


3-7 


with 


iVi(7) = 


47r 


iV2(7) = 


Nai< a) = N 2 {'y) nor^ 

22-77r3/2 p,'^ 


3-7 


(B8) 


(B9) 


3-7’ 3-7 r(7-l/2)- 

Values of and N 2 are given in the table. For 7 > 3/2 the two quantities are very similar; they begin to depart 


7 

1 

3/2 

7/4 

2 

5/2 

Ni 

6.28 

8.38 

10.05 

12.56 

25.13 

N2 

3.14 

6.98 

9.40 

12.56 

26.17 

N 2 /N 1 

0.50 

0.83 

0.94 

1.00 

1.04 


significantly for smaller 7. 


APPENDIX 

Anomalous Relaxation 


A Monte- Carlo algorithm for d escribing the evolution of L in the “anomalous relaxation” (L < Lsb) regime was 
presented in iMerritt et al.l (j2011h . That algorithm was based on a simple model for the time- and space-dependence 
of the '/N perturbing potential. Following is an analytic derivation of the angular momentum transition probabilities, 
and the corresponding diffusion coefficients, that are implied by the same Ha miltonian model. The results of the 
deriva tion presented in this appendix were the basis for the functional forms that iHamers. Portesfies Zwart fc Mer~rittl 
(|2014f) fit to diffusion coefficients extracted from their A^-body data, and which appear here in il3.3l 
Consider a star orbiting in the potential 

$(r) = $Kepler + ‘I’n + ■ (Cl) 

Here, <i>Kepier = —GM./r is the (Newtonian) potential due to the SBH; <i>N is the potential from the spherically- 
distributed mass; and $^ is the potential due to the y/N asymmetries in the stellar distribution. We assume that 

I'^nI “C |$Kepler|- 

If the mass density falls off as a power of radius, p{r) = po{r/ro)~'^, then 

= (2 T) (;/) + 

for 7 ^ 2; for 7 = 2 the dependence of <I>n on radius becomes logarithmic. Following IMerritt et al.l (120111) . we assume 
that the perturbing potential, as experienced by a test star of semimajor axis o, is given by 

=—aS{a) cos6 =—aS{a)-, S{a) = — . (C3) 

We are assuming for the moment that $^7^ is independent of time. Expressing the two perturbing potentials in 
Delaunay variables and averaging over the unperturbed (Keplerian) motion yields 


- GM*(r < a) 
<Pn = 


(2 - 7)0 
$^77 = aS'(a)e sin f sin w. 


(1 -I- ai - q; 2 ^^) 


(C4a) 

(C4b) 


Here, i = Vl — e^, i = 7r/2 corresponds to the x — z plane and w = 7r/2 describes an orbit that is elongated along z. 
The e xpression for $n assumes £ <C 1; the quantities ai(7), a2(7) are both of order unity and are given in iMerrittI 
(IMl §4.4.1). 

Ignoring constant terms (including terms that depend only on a), the averaged Hamiltonian is then 


H = -5- = —£ ^ — An£^ + a e sin i sin w. 


An = 


a2 


M^{a) 


3(2-7) M. Tg 


1 S {a) a 

2GM,la?Vg 


(C5a) 


(C5b) 
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and 




9 ’ 


-P(fl) a 


Discarding terms of order or smaller, the averaged Hamiltonian becomes 


(C6) 


(C7) 


H =—£ +H^sinisinw. 

Note that the “mass precession” terms vanishes to this order in i. The first term in equation (IC7I) represents GR 
precession; the second represents the effects of the VN torques. 

Writing r = vot, Hamilton’s equations of motion for the osculating elements are 

M _ dH 
dr doj 

(C8) 


dw dH 
dr dr 
dfl dH 

dr d£. 


= £-\ 
= -A 


sini ’ dr 

Henceforth setting sini = 1, the dependence of £ on a; is given by 

1 


= — A sini cos w, 
= 0 . 


[(£2 - £ 1 ) sinw + (£1 + £ 2 )] 


2£i£2 


(C9a) 

(C9b) 


where {£i,£ 2 } are the extreme values of £ and go = —H/A^ is the “energy.” The equation of motion for £ can be 
expressed in terms of £ alone as 


M 

dr 


= 


1 1 

Ti~J 


1 1 

Ti~J 


with a similar expression for uj, and the full precessional period is 

27r go 


To = 


(sS - 


where 


thus go = (£2 +£i )/(£2 -£i) and 


(50 + 1) \ 4 — A }- {go - 1) 


-1 


Vn 
£1 + £2 


— £a.v 


An 
1 go 

Ay/N 5o “ 1 


(CIO) 

(Cll) 

(C12) 

(C13) 


Solutions obtained so far describe “coherent resonant relaxation:” changes in a star’s angular momentum for times 
shorter than the coherence time. Now, suppose that the orientation of the torquing potential changes, instantaneously, 
at random times separated by £coh- Of course this is a crude oversimplification since in reality the torquing potential 
is changing gradually; on the other hand, for stars near the Schwarzschild barrier, the GR precession time is expected 
to be comparable to the coherence time (equation [BS]). 

Let the angle between the new z-axis , and the orbital semimajor axis at the moment of the switch, be loi. At this 
moment, £ has the value is, and go changes to gi, where 


9i — ^/{A^is) — sinwi. 

Since the probability distributions of r and wi are uniform, we can write 

P {is,gi) disdgi = P (r, wi) drduji = 


27rTn 


(014) 


(CIS) 


The factor four on the RHS accounts for the fact that £ varies over its full range in one-fourth of a precessional 
period. Then 


^ 2 5(t,wi) _ 2 ^ 

P\^Si9l) rp f\( p \ - T' 

TTTod{ts,gi) ttTo 
dis dgi 


j-i^d{is,gi) 


d{T,U}l) 


dr dtoi 


J = A^\J 1 - [go - 1/ {A^is)] \Jl-[gi- l/(A^£s)] . 


(C16a) 

(C16b) 

(C16c) 


















25 


9o = 



Fig. 10.— P{gi\go) computed in two ways, for go = 3. 


Integrating over is gives the probability distribution for gi. This can be written 


-P(5i|5o) 


dx 




x ^ yj(x — Xi ){ x 2 — x )[2 — {x — Xi)] [2 — ( x 2 — x)] 


(C17) 


where 

for 5i < 5o, and 


xi = go - X2 = gi + 1 

xi = gi - 1 , X2 = go + l 


(C18) 

(C19) 


for gi > go- Figure [TO] shows a numerical integration of equation (jC17ll (solid line) for go = 3, compared with the 
results of Monte-Carlo experiments based on the equations of motion (|C8I) . We note here the asymmetry of the derived 
transition probability. 

The diffusion coefficients are 


1 /-So+S 

{{Ago)^) = -— / {gi-go) P{gi\go)dgi. (C20) 

-^coh JgQ — 2 

The integral (IC20I) can be broken into two pieces, /> and /<, corresponding to go > gi and go < gi respectively. In 
the case of {Ago), 


pXi-\- 2 px 

(I'To^^T'cohj ^>= J ( 2^2 - ail - 2) dx2 J 


dx 


\J { X 2 - x){x - Xi)[2 - { X 2 - a;)] [2 - {x - a;i)] 


- 2 


(1 — w)dg 


5o “ 1 Jo \/ w {2 - w ) {w + go - I) 


■ sin 


-1 


2 — w 


and a similar calculation gives 


( 2 /< - 2 _ ^ + 2 ^ 


(I — w)dw 


5o “ 1 Jo \/ w {2 - w ) {-w -1-50 + 1) 


• -1 i^-w 
■ sm ■ ' 


The sum is 


(/< +/>) = -^/(5o), 
^ ^ ^ 9q 


(5q — a:2)2\/l — x'^ 


sm 


I - a; TT^ 5 g 


2 Hgl-lf'^' 


(C21a) 

(C2Ib) 

(C22) 

(C23a) 

(C23b) 
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Thus 


(Ago) = - 


5o TrTcohTb^^ 


Vn 


-I{go) = -T^^- 

^coh 90 


(C24) 


Proceeding as before to evaluate ((Ago)^): 

((Ago)^) = 

K{9o) = 


irToA^^Tcoh 


Vn 

5 , (ag - 1)“-'^ 

2 go 


{!< + I>) — If ,— K{go), 

^ coh 


<?o 


(C25) 


The function K{go) is almost independent of go: 

X(3) = 1.042, iC(5) = 1.015, K{8) = 1.006, ^(10) = 1.004, K{20) = 1.001, K{oo) 

so that 


= 1 


1 


To a good approximation then, 

and the time scales for changes in go are 


((Ago)^) = (1.0™)— . 

^coh 


(Ago) = ((Ago)^) = ^ 

^coh 90 ^coh 


9o^coh- 


(Ago) 

-1 

((Ago)^) 

go 


5o 


Identifying £ with £av (equation IC13P , diffusion coefficients in £ become 

3,. 


{A£)^^ 

T 


1 + 2 

((A£)^)«— l + 2[A^ £) 


(C26) 

(C27) 

(C28) 

(C29a) 

(C29b) 


where T(a) = tcoh{a)/. 
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